home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / ZBRENT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  87 lines
  1. FUNCTION zbrent(x1,x2,tol: real): real;
  2. (* Programs using routine ZBRENT must externally define a function
  3. fx(x:real):real a root of which is to be found. *)
  4. LABEL 99;
  5. CONST
  6.    itmax=100;
  7.    eps=3.0e-8;
  8. VAR
  9.    a,b,c,d,e: real;
  10.    min1,min2,min: real;
  11.    fa,fb,fc,p,q,r: real;
  12.    s,tol1,xm: real;
  13.    iter: integer;
  14. BEGIN
  15.    a := x1;
  16.    b := x2;
  17.    fa := fx(a);
  18.    fb := fx(b);
  19.    IF (fb*fa > 0.0) THEN BEGIN
  20.       writeln('pause in routine ZBRENT');
  21.       writeln('root must be bracketed'); readln
  22.    END;
  23.    fc := fb;
  24.    FOR iter := 1 TO itmax DO BEGIN
  25.       IF (fb*fc > 0.0) THEN BEGIN
  26.          c := a;
  27.          fc := fa;
  28.          d := b-a;
  29.          e := d
  30.       END;
  31.       IF (abs(fc) < abs(fb)) THEN BEGIN
  32.          a := b;
  33.          b := c;
  34.          c := a;
  35.          fa := fb;
  36.          fb := fc;
  37.          fc := fa
  38.       END;
  39.       tol1 := 2.0*eps*abs(b)+0.5*tol;
  40.       xm := 0.5*(c-b);
  41.       IF ((abs(xm) <= tol1) OR (fb = 0.0)) THEN BEGIN
  42.          zbrent := b; GOTO 99 END;
  43.       IF ((abs(e) >= tol1) AND (abs(fa) > abs(fb))) THEN BEGIN
  44.          s := fb/fa;
  45.          IF (a = c)  THEN BEGIN
  46.             p := 2.0*xm*s;
  47.             q := 1.0-s
  48.          END ELSE BEGIN
  49.             q := fa/fc;
  50.             r := fb/fc;
  51.             p := s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
  52.             q := (q-1.0)*(r-1.0)*(s-1.0)
  53.          END;
  54.          IF (p > 0.0) THEN  q := -q;
  55.          p := abs(p);
  56.          min1 := 3.0*xm*q-abs(tol1*q);
  57.          min2 := abs(e*q);
  58.          IF (min1 < min2) THEN min := min1 ELSE min := min2;
  59.          IF (2.0*p < min) THEN BEGIN
  60.             e := d;
  61.             d := p/q
  62.          END ELSE BEGIN
  63.             d := xm;
  64.             e := d
  65.          END
  66.       END ELSE BEGIN
  67.          d := xm;
  68.          e := d
  69.       END;
  70.       a := b;
  71.       fa := fb;
  72.       IF (abs(d) > tol1) THEN BEGIN
  73.          b := b+d
  74.       END ELSE BEGIN
  75.          IF (xm > 0) THEN BEGIN
  76.             b := b+abs(tol1)
  77.          END ELSE BEGIN
  78.             b := b-abs(tol1)
  79.          END
  80.       END;
  81.       fb := fx(b)
  82.    END;
  83.    writeln('pause in routine ZBRENT');
  84.    writeln('maximum number of iterations exceeded'); readln;
  85.    zbrent := b;
  86. 99:   END;
  87.